Critical behaviors of sheared frictionless granular materials near jamming transition 
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Critical behaviors of sheared dense and frictionless granular materials in the vicinity of the jam- 
ming transition are numerically investigated. From the extensive molecular dynamics simulation, 
we verify the validity of the scaling theory near the jamming transition proposed by Otsuki and 
Hayakawa (Prog. Theor. Phys., 121, 647 (2009)). We also clarify the critical behaviors of the shear 
viscosity and the pair correlation function based on both a mean field theory and the simulation. 
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I. INTRODUCTION 

Let us consider mechanical properties of grains which 
are packed in a container. When the density is low 
enough and the effect of gravity is negligible, there is no 
pressure acting on the wall of the container. However, 
when the density exceeds a critical value, the pressure 
acting on the wall becomes finite. Such kind of transi- 
tion for the rigidity or the stress is known as the jamming 
transition. 

Jamming is a key concept to characterize the transi- 
tion of athermal systems such as granular materials ^1], 
colloidal suspensions emulsions and forms 0. Liu 
and Nagel suggested that a unifying description might 
be possible to cover both glass transition in thermal sys- 
tems and athermal jamming transition Indeed, there 
are many similarities between the jamming and the glass 
transition. For exaniple, the rheological properties of the 
jammed systems fs', ^ are similar to those of glassy ma- 
terials and the granular materials near the jamming 
transition point (point J) exhibit the dynamical hetero- 
geneity which is one of the characteristics of glassy mate- 
rials [6, 8; 9, l3, EH, 13]. From these similarities, the 
mode coupling theory for glassy materials [3, [l^, [3| is 
applied to granular materials, but it fails to describe the 
jamming transition [iTj . 

Recent studies have revealed that the jamming is the 
transition to appear in the pressure, the coordination 
number, the elastic moduh [IS, |19|], and the soft modes 
[2^ . In particular, Olsson and Teitel 0], and Hatano [2lj 
numerically found scaling relations for the shear stress, 
the pressure, and the kinetic temperature, which are 
characterized by some critical exponents. A similar scal- 
ing relation is observed in the simulation of Josephson 
junction arrays [l^. These scaling relations indicate that 
jamming is a continuous transition in which the stress is 
zero at the critical point. In order to understand the 
critical behaviors of the jamming transition, it is impor- 
tant to determine the critical exponents. In the previous 
paper [2^, the authors theoretically obtained the criti- 
cal exponents for the scaling relations of the shear stress, 
the kinetic temperature and the characteristic frequency. 



which characterizes the collisional energy loss, for sheared 
frictionless granular materials. They also numerically 
verified the validity of their theoretical predictions for 
the linear spring model 23]. 

In this paper, we will verify the validity of the predic- 
tions in Ref. [23| for sheared frictionless granular ma- 
terials in various situations. We will also discuss the 
behaviors of the viscosity and the pair correlation func- 
tion in details. In the next section, we will summarize 
the theoretical predictions by the present authors [23| . 
In Sec. IIIH we will compare the results of our exten- 
sive simulations with the theoretical predictions. Section 
mil consists of five subsections: Section fill Al will be de- 
voted to the explanation of the setup of our simulations, 
in Sec lIIIB] we will present various scaling plots to ver- 
ify the scaling theory, in Sec. IIII CI we will discuss the 
force law dependence of the critical exponents, we will ex- 
plain the density dependence of critical variables in Sec. 
IIIIDI We will explain the results for nearly elastic cases 
in Sec. IIII El Section HVl will be devoted to the explana- 
tion of critical properties of the pair correlation function 
in the vicinity of the jamming transition, which was not 
discussed in Ref. |23|. In Sec. El we will discuss and 
conclude our results. In Appendix A, we will summarize 
the method for the theoretical prediction which contains 
some generalized results beyond Ref. ^23'] based on a dif- 
ferent method for the derivation. In Appendix B, we will 
derive some relations, which are necessary to discuss the 
critical behavior of the pair correlation function. 



II. MEAN FIELD THEORY 

In this section, we briefly summarize the theoretical 
results in Ref. [2^ . Let us consider D-dimensional gran- 
ular assemblies under an uniform shear with shear rate 
7. The system includes N grains, each of which has the 
identical mass m. The packing fraction of the system 
and the critical fraction at point J are respectively de- 
noted by (j) and 0j. Throughout this paper, we assume 
that granular particles are frictionless, where any contact 
force acts along the line to connect two centers of mass of 
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contacting grains. In most of our arguments, we assume 
that the elastic interaction between the grain i located 
at ri and the grain j at Tj is given by 



fc\{rij) = fee [aij - rij) {(Jij - rij) 



A 



(1) 



where k and the elastic constant and the distance 

between the grains rij = \rij\ = \ri — rj\, respectively. 
aij = {ai + aj)/2 is the average of the diameter of the 
grain i and the grain j with diameters <Ti and aj. Q{x) is 
the Heaviside step function satisfying 8(x) — 1 for x > 
and <d{x) = for otherwise. Many researchers use the 
linear spring model (A = 1) associated with the viscous 
contact force 



(2) 



where 77 is the viscous parameter. Here, Vij^n is the 
relative normal velocities between the grains given by 
Vij^a = (vi — Vj) ■ Tij/rij, where Vi and Vj are the veloc- 
ities of the grain i and the grain j, respectively. On the 
other hand, the Herzian model (Eq.([T]) with A = 3/2) 
associated with the corresponding viscous contact force 

I vis{f ij ~ —7/0 {(Tij — rij) (Jij ~ TijVij^ji (3) 

is more realistic one for three-dimensional spheres. 

In the previous paper [23|, we introduced the scaling 
relations 
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(4) 
(5) 
(6) 
(7) 



where $ = — (pj. The kinetic temperature T, the shear 
stress S, and the pressure P are respectively given by (23 | 
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where we have introduced Pj = m[vi — c(rj)] with the 
average velocity Ca{r) = jySa^x at the position r, V is 



the volume of the system, and (•) represents the ensemble 
average. The characteristic frequency lo is defined by 



(11) 



where n is the number density. This to is reduced to the 
collision frequency in the unjammed phase determined 
by the balance between the viscous heating and the col- 
lisional energy loss. T^{x), 5+ (a;), V+{x) and W+ix) 
are the scaling functions above point J (in the jammed 
phase). On the other hand, the scaling functions below 
point J (in the unjammed phase) are given by T-{x), 
S-{x), V-{x) and >V-(x). Equations ©-(III) contain 
the critical exponents a:^, x^, y^, y^, y'^, y!y, z^, and z^. 
It should be noted that At^d, Aw,d, ^p,_d, tu, 

Sd, wd, and pd are the constants depending only on the 
dimension D. 

From the scaling theory explained in Appendix A, 
there are the scaling relations in the unjammed phase 
in the limit 7 — > as 

T ^2|(j,|a:^(l-2/a:T)^ 5" ^ ^2 |(j,|a^(l-2/yT) ^ 

p ^ j^\<^\v'4^-^/yy\ a; - 7|$|^*(i-i/^^). (12) 
Similarly, T, 5, P, and lo satisfy 

P ^ uj - |$|^*, (13) 

in the zero shear limit of the jammed phase. The scaling 
relations at point J, i.e. (f> ~ 0, are given by 



p. 



7 . 



(14) 



As explained in Appendix A and Ref. [23| , the critical 
exponents are given by 

a A 2A + 4 ^ 2A 

X^ = 2 + A, 2;^ = ^—-, y^=A, y^ = -^—^ 

/ A , 2A A A , , 

y^ = A, y^^ z$ = -, z., = (15) 



A + 4' 



Note that the scaling of the pressure and the determina- 
tion of the exponents y'^ were not discussed in Ref. [1^ . 
The derivations which are a little different from the orig- 
inal one are explained in Appendix A. We should stress 
that the exponents are independent of the spatial dimen- 
sion D, and the characteristic feature of the jamming 
transition is the A-dependence of the critical exponents. 
Indeed, O'Hern et al. found that the pressure behaves as 
P ^ |$|^ for unsheared jammed systems in the vicinity 
of the jamming point. Hatano also indicated that the 
critical exponents for the sheared granular materials dif- 
fer in the cases of A = 1 and 3/2. Explicit A-dependence 
of the exponents in Eq. (fTS)) except for y'^ for the sheared 
granular materials was obtained in Ref. [1^. It should 
be noted that our theory can be generalized to the case 
of the contact force 



/ol(rij) = e(cry - r,j)T{nj). 



(16) 
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Here, the exponent A in Eq. ([15)) should be determined 
from the relation \iinr--^a-j T{rij) oc (cr.y — Vij)^ for this 
case. If we use an analytic function T{rij) such as the 
repulsive Lennard- Jones force 
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(17) 



because we can use the 
-.,)+^"(fT,,)/2(a,, - 



the value of A should be A = 1 
expansion T{rij) ~ T' {aij){aij — 

)^ + ' ' ■ ■ Note that Eq. p7)) differs from the usual 
form for the Lennard- Jones potential in order to satisfy 



III. COMPARISON BETWEEN THEORY AND 
SIMULATION 

In this section, let us verify the validity of the theory 
[2^ 1 by the molecular dynamics simulation. As stated in 
Introduction, we have already confirmed that all of our 
results of the simulation for A = 1 under one size distri- 
bution of grains with fixing the strength of inelasticity are 
consistent with the theory [235 , but nobody has checked 
its validity for general A under various situations. Thus, 
we will perform extensive simulations with changing A, 
the size distribution, and the inelasticity etc. 

In the first part pil Ap . we explain the details of our 
model. In the second part piIBp . we demonstrate the 
validity of the scaling theory given by Eqs. (H])-© with 
Eq. (|15p from the test of scaling plots under various 
conditions. In the third part (jlll C|) . we check the the- 
ory for the A-dependence of the critical exponents. In 
the fourth part piIPp . we discuss ^-dependences of the 
characteristic frequency oj in the jammed region and the 
viscosity ^ = S/j in the unjammed region. In the last 
part pil E[) . we explain the behaviors for nearly elastic 
cases, while the results in the first five parts are obtained 
for strongly dissipative cases. 



call the monodisperse system consists of only one type of 
particles, whose diameters are ctq. The critical density 
of the monodisperse system is believed to be 0.639 for 
D — 3 from a numerical simulation [l^ . 

The time evolution equations of the position and 
the velocity Vi of the ith particle are given by 



drj 

dt 

dv, 
I—— 
dt 



(18) 



X! {/oi(''y ) + /vis(nj, %■,„)} — , (19) 



where the elastic force fci{rij) is given by Eq. ([T]) for 
most of cases except for the case of the repulsive Lennard- 
Jones potential. fvis{rij,Vij,n) given by Eqs. ([2]) or ([3]) is 
the viscous contact force between particles. In order to 
realize an uniform velocity gradient 7 in y direction and 
macroscopic velocity only in the x direction, we adopt 
the Lees-Edwards boundary conditions. 

In our simulation with the elastic contact force given by 
Eq. Il]), TO, (To and k are set to be unity, and all quantities 
are converted to dimensionless forms, where the unit of 

time scale is mal^'^/k. In the case of the contact force 
given by Eqs. ^ and e in Eq. dTT]) is set to be 

unity instead of k, and the unit of time scale is ma'^/e. 
We adopt the elastic constant fc = 1.0 or e = 1.0, and 
the viscous constant 77 = 1.0 for most of cases except for 
nearly elastic cases in Sec. IIII El This situation corre- 
sponds to the constant restitution coefficient e = 0.043 
for the linear spring model. We use the leap-frog algo- 
rithm, which is second-order accurate in time, by using 
the time interval At = 0.2 for the cases of the contact 
force given by Eq. ([1]) with checking the convergence un- 
til At = 0.05. In the simulation with the contact force 
given by Eqs. ^ and dTT]), we use M = 0.01. 



B. Scaling plots 



A. Setup 

We examine three different systems on dispersion of 
diameters of grains. The first system we call the poly- 
disperse system consists of four types of grains, and the 
diameters of grains are 0.7cto, O.Stroj 0.9ao and ctoj where 
the number of each type of grains is iV/4. The poly dis- 
perse system has been studied in Ref. [235, where the 
critical fraction is estimated as = 0.84285 for D = 2, 
= 0.64455 for D = 3 or = 0.4615 for Z) = 4 [H. 
The second system which we call the bidisperse system 
consists of two types of grains. The diameters of grains 
are 5ao/7 and ctq, where the number of each type of grains 
is N/2. The bidisperse system has been studied by many 
researchers [H, IH, [2^ , and the critical density j 
is known as 0.648 for D = 3 by the numerical simula- 
tion of static granular packings [Iq . The final system we 



Figures [T]|3] show the scaling plots of the polydisperse 
system with A = 3/2 based on Eqs. with Eq. 

for the various dimensions D = 2,3, and 4. (See 
Ref. [23\ for the scaling plots of the polydisperse system 
with A = 1). We use Eq. ^ for the viscous contact 
force. Here, the number of the particles N is 2000, and 
the shear rate 7 is ranged between 5 x 10"'' and 5 x 10~^ 
for D = 2,3 and between 5 x 10^^ and 5 x 10^"* for 
D — 4. We also use the amplitudes and the adjustable 
parameters {to, At.o, sd, As,d,Pd, Ap^D,WD, A^^d) = 
(0.01,15.0,0.03,0.02,0.02,0.3,0.1,0.12) for D = 2, 
(0.01,6.0,0.04,0.03,0.025,0.3,0.15,0.25) for D = 3, 
(0.1, 0.45, 0.05, 0.05, 0.04, 0.5, 0.1, 0.45) for £> = 4. AU of 
the data converge to the universal master curves. These 
results verify the validity of our theoretical predictions in 
Eq. (Uni). 

In Figs. [H] and [SI we show the scaling plots of S for 
the bidisperse system and the monodisperse system with 
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FIG. 1: (Color online) Collapsed data of the shear rate depen- 
dence of the kinetic temperature T in the polydisperse system 
(N = 2000) with A = 3/2 using the scaling law, Eq. (gj, 
for 73 = 2, 3 and 4. The dashed line, the dotted line and the 
solid line are proportional to 7, 7^ and 7^^ , respectively. The 
legends show the dimension D and the volume fraction (j) as 
{D,(j)). The critical exponents estimated from Eq. (|15fl are 
= 9/2 and Xj = 14/11. 



FIG. 3: (Color online) Collapsed data of the shear rate 
dependence of the pressure P for the polydisperse system 
(iV = 2000) with A = 3/2 using the scaling law, Eq. ©, 
for D = 2, 3 and 4. The dotted line and the solid line are 
proportional to 7^ and 7^t, respectively. The legends show 
the dimension D and the volume fraction </!> as {D,(j>)- The 
critical exponents estimated from Eq. (|15[) are y'^ — 3/2 and 

= 6/11. 
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FIG. 2: (Color online) Collapsed data of the shear rate de- 
pendence of the shear stress S for the polydisperse system 
{N = 2000) with A = 3/2 using the scaling law, Eq. ©, 
for D = 2, 3 and 4. The dotted line and the solid line are 
proportional to 7^ and 7^^, respectively. The legends show 
the dimension D and the volume fraction (j> as {D,(j>). The 
critical exponents estimated from Eq. (I15|l are — 3/2 and 
y-r = 6/11. 



FIG. 4: (Color online) Collapsed data of the shear rate de- 
pendence of the characteristic frequency co for the polydisperse 
system (A'^ = 2000) with A = 3/2 using the scaling law, Eq. 
(|7|), for D = 2,3 and 4. The dotted line and the solid line 
are proportional to 7 and 7^^ , respectively. The legends show 
the dimension D and the volume fraction (f> as {D,(j}). The 
critical exponents estimated from Eq. (|15|l are = 3/4 and 
= 3/11. 



A = 1 and D = 3, respectively, where both systems in- 
clude N — 4000 particles. The viscous contact force is 
given by Eq. and the shear rate 7 is ranged be- 

tween 5 X 10^^ and 5 x 10^'"'. The amplitude and the 
adjustable parameter are {sd,As,d) = (0.035,0.04) for 
the bidisperse system and (0.025, 0.035) for the monodis- 
perse system. These scaling plots support the validity 
of our prediction. The scaling plots for T, P, and lu for 
the bidisperse and the monodisperse systems also exhibit 
elegant scalings, but we omit these figures in this paper. 

We have checked the validity of our scaling theory in 
larger systems. Figure [7] shows the scaling plot of S in the 
three-dimensional monodisperse system with N = 20000 
particles, where the shear rate 7 is ranged between 5 x 
10~^ and 5 x 10~^. The parameters and the guide lines 
are the same as those for Fig. [51 This scaling supports 
the validity of our theory even in the larger system. 
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FIG. 5: (Color online) Collapsed data of the shear rate 
dependence of the shear stress S for the bidisperse system 
(A'^ = 4000) with A = 1 using the scaling law, Eq. ([5}, for 
D = 3. The dotted line and the solid line are proportional to 
7^ and 7''"' , respectively. The legends represent the volume 
fraction (p. The critical exponents estimated from Eq. psp 
are y^ = 1 and y^ = 2/5. 
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FIG. 6: (Color online) Collapsed data of the shear rate de- 
pendence of the shear stress S for the monodisperse system 
(A'^ = 4000) with A = 1 using the scaling law, Eq. for 
D = 3. The dotted line and the solid line are proportional 
to 7^ and 7^^, respectively. The legends show the volume 
fraction (p. The critical exponents estimated from Eq. (|15[) 
are = 1 and y-, = 2/5. 



FIG. 8: (Color online) Collapsed data of the shear rate de- 
pendence of the shear stress S in the monodisperse system 
(A'^ = 2000) with the contact force given by Eqs. (|16p and 
(fT7|) using the scahng law, Eq. ([5} for 7? = 3. The dotted 
line and the solid line are proportional to 7^ and 7^^" , respec- 
tively. The legends show the volume fraction (j>- The critical 
exponents estimated from Eq. (|15|) with A = 1 are y^ — 1 
and y-y — 2/5. 
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FIG. 7: (Color online) Collapsed data of the shear rate depen- 
dence of the shear stress S in the larger monodisperse system 
(iV = 20000) with A = 1 using the scaling law, Eq. for 
D = 3. The dotted line and the solid line are proportional 
to 7^ and 7^^, respectively. The legends show the volume 
fraction (p. The critical exponents estimated from Eq. (|15|) 
are 2/0 = 1 and y^ = 2/5. 



In order to verify the validity of our theory in more 
general cases than that of Eq. ([1]) , we examine the scaling 
plot for 5* in the three-dimensional monodisperse system 
with the elastic contact force given by Eqs. and (fT7|) 
and the viscous contact force given by Eq. ^ in Fig. 
[HI Here we use the number of the particles N = 2000, 
and the shear rate 7 is ranged between 10~^ and 10~^. 
The amplitude and the adjustable parameter are given by 
{sd,As^d) — (0.004,3.0), where we adopt the exponents 
given by Eq. with A = 1. The scahng in Fig. [5] 

supports the validity of our theory for the elastic contact 
force given by Eq. ITF 



C. A-dependence of the critical exponents 

In order to verify the A-dependence of the critical ex- 
ponents, we plot y-y and y'^ versus A for the polydisperse 
system with D = 2, N = 4000, and the viscous force 




FIG. 9: Plot of y-y versus A for the polydisperse system (A'^ : 
8000) with D = 2. 



given by Eq. ^ in Figs. [51and[TUl respectively. Here, 
in order to obtain Figs. [9land[T0l we have estimated y-y 
and y^j,, respectively, from the shear stress <S'(7, 0) and 
the pressure P{j, (j)) at (j) = (j),/ as 



log(^(7i>j))-log(g(72,0j)) 
log7i-log72 
, ^ log(P(7i,0j))-log(P(72,0j)) 
^ Iog7i-log72 



(20) 
(21) 



for (71,72) = (5 X 10-^1.5 X 10-6)^(1.5 x 10^6,5 x 
10-6), (5 X 10-6, 1.5 X 10-5), (1.5 X 10-^ 5 X IQ-^), (5 x 
10-^5 X 10-6), (5 X 10-6,5 X 10-5), and plotted the av- 
erages of y-y and y'^ for different (71,72) with the error 
bars, whose lengths are twice of the standard deviation 
of y-y and y',y. The exponent y^ reasonably agrees with 
the prediction Eq. (fTS]) in the wide range of A. Although 
y'-y is a little deviated from the theoretical prediction, we 
believe that the deviation becomes smaller if we use the 
data for smaller 7 and larger TV. 

Figures [TT] and [T^ demonstrate whether the exponents 
for 5' predicted by Eq. can be used for the scaling 
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FIG. 10: Plot of y'^ versus A for the polydisperse system 
{N = 8000) with D = 2. 
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FIG. 11: (Color online) CoUapsed data of the shear rate 
dependence of the shear stress S for the polydisperse system 
(iV = 8000) with A = 0.5 using the scaling law, Eq. ©, for 
D — 2. The dotted line and the solid line are proportional 
to 7^ and 7^^, respectively. The legends show the volume 
fraction (p. The critical exponents are estimated from Eq. 
(fTS]) as j/0 = 1/2 and y^ = 2/9. 
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FIG. 12: (Color online) Collapsed data of the shear rate 
dependence of the shear stress S for the polydisperse system 
(iV = 8000) with A = 2.0 using the scaling law, Eq. ©, for 
D = 2. The dotted line and the solid line are proportional 
to 7^ and 7^T, respectively. The legends show the volume 
fraction (j>- The critical exponents are estimated from Eq. 
psp as j/0 = 2 and = 2/3. 
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FIG. 13: Plots of u! versus $ of the polydisperse systems for 
D = 2 and A = 3/2 with 7 = 5 x 10"®, 5 x 10"'', and 5 x 10"*. 



plots of A = 0.5 and 2, respectively. The shear rate 7 is 
ranged between 5 x 10~^ and 5 x 10~^. The amplitude and 
the adjustable parameter are (sD,^s,n) — (0.015,0.05) 
for A = 0.5 and (0.05,0.02) for A = 2.0, respectively. 
These scalings in Figs. [11] and [T^] as well as the evalu- 
ated exponents by Eqs. (^0]) and (pij) strongly support 
the validity of the theoretical predictions in Eq. for 
arbitrary A. Thus, our theory can be used for any A. 



D. ^-dependence of critical variables 

Here, we examine ^-dependences of quantities pre- 
dicted in Eqs. ^ and ^ with Eq. First, we plot 
uj versus |$| in the jammed phase for the two-dimensional 
polydisperse systems {N ~ 2000) with A = 3/2 in Fig. 
[131 We adopt Eq. ^ for the viscous contact force. 
Note that the corresponding results for A = 1 have been 
reported in [2^. As the shear rate 7 decreases, uj ap- 
proaches CO ~ as predicted in Eq. (|13p with Eq. 
P^ . There is a plateau in the small 7 region, but the 
value of it decreases as 7^^ in the limit 7^0, which can 



be predicted from Eq. ([7]) because yV+{x) ~ x^"/ with 
X = 7/|(I)|^<<>/^T- 00. 

Next, we show ^-dependence of the viscosity /i = 5/7 
in the unjammed phase, which is predicted as /i ^ 
7|$|-4 ^ |$|-4 gqg^ ^ g^j^j ^ ^Q^g ^Yi&t 

the critical exponent for fi is independent of A and D. 
Figure [HI includes the data of ^/j as a function of |$| 
for the polydisperse system {N = 2000) with A = 1 and 
3/2. "We respectively adopt Eqs. ([2]) and ([3]) for the vis- 
cous contact forces in in the cases of A = 1 and 3/2. 
Both of the data for A = 1 and 3/2 satisfy the theoret- 
ical prediction /i/7 ~ as predicted, although there 
is a plateau when |$| 0. The value of the plateau 
proportional to 7^^"^ can be understood by Eq. ([5|) 23|. 

There are some previous studies on the diver gen ce of 
the viscosity /i. For example, Losert et al. [23] ob- 
served exponents larger than 1 from their experiment. 
The critical exponent of the shear viscosity for foams 
and emulsions in Ref. Q is between 1 and 2. For col- 
loidal suspensions, the viscosity is believed to diverge 
as |29|. Garcia-Rojo et al. [H] reported that 

the scaled viscosity diverges as i^/VT ^ {(f)^ — 4>)~'^ for 



FIG. 14: Plots of versus <1> of the two-dimensional 
polydisperse systems for and A = 1 and 3/2 with -y — 
5 X 10"^5 X 10"*', and 5 x 10"^ 



FIG. 16: The snapshot of the two-dimensional monodisperse 
system {N = 2401) with A = 1, it = 1.0, and rj = 0.00225. 
We use Eq. ([2]) for the viscous contact force. The shear rate 
is 7 = 5 X 10~^, and the density is ^ = 0.80. 
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FIG. 15: m/VT as a function of <!> for the three-dimensional 
btdisperse system (A'^ = 4000) with A = 1 and 7 = 5 x lO"'^, 
where the solid line and the dashed line are proportional to 
{4>fj, — (f))'^ and {4>j — (f>)~^, respectively. 



(j)^ < (j)j in two-dimensional elastic monodispersed hard- 
disks. This result is contrast to our prediction that the 
viscosity scaled by the temperature ^/T diverges at point 
J as /i/\/T ^ {(j),] — cj))~'^ , obtained from Eq. ([T2l) with 
Eq. (HID. 

In order to clarify whether our prediction is valid for 
sheared frictionless granular materials in the vicinity of 
the jamming transition, we examine the possibility that 
the viscosity satisfies n/y/T ~ (0^ — 0)^-^ with a fitting 
parameter (f)^ [28^ in Fig. [15] for the bidisperse system 
{N = 4000) with D = 3, A = 1 and 7 = 5 x 10~^. 
Actually we can fit the data of our three-dimensional 
simulation by I-l/Vt ~ (0^ — 0)^"'^, but the viscosity is 
still finite even for </) > 0^ = 0.632. We, thus, conclude 
that the viscosity /i in the sheared granular materials 
does not satisfy fi/ VT ^ (0^ — and the behavior of 

/i is consistent with our prediction. Here, we should note 
that our theory is no longer valid for the two-dimensional 
monodisperse system, where the shear band occurs (Fig. 
I16p . which is not assumed in our theory. 

Berthier and Witten [2^, |3^ reported that the relax- 
ation time in the zero-temperature limit of the three- 



dimensional equilibrium bidisperse system diverges at 
(j)G — 0.635, which is smaller than — 0.639. Although 
this idea might be attractive to characterize universal 
feature of the jamming transition, we could not find such 
divergence, as shown in Fig. 1151 for frictionless sheared 
granular materials near the jamming transition. 



E. Nearly elastic cases 

One might think that our scaling theory is only valid 
when the dissipation of the system is strong enough. In- 
deed, we have used the viscous constant rj — 1.0, which 
corresponds to the restitution coefficient e — 0.043 for 
the linear spring model. In order to check the validity of 
our theory in the nearly elastic system (e ~ 1), we per- 
form the numerical simulation of the three-dimensional 
monodisperse system with A = 1 and rj = 0.018, which 
corresponds to the restitution coefhcient e = 0.96. We 
use the viscous contact force in Eq. ^ with the number 
of particles N — 2000. The shear rate 7 is ranged be- 
tween 5 X 10^^ and 5 x 10^^. The amplitude and the ad- 
justable parameter are given by {sd,As,d) — (0.03, 0.04). 
The scaling plot of S in this system is shown in Fig. 1171 
This figure supports the validity of our scaling even in 
the nearly elastic system. 



IV. PAIR CORRELATION FUNCTION 

In this section, we discuss the behaviors of the spatial 
correlations in the vicinity of the jamming transition. In 
particular, we focus on the critical behaviors at the first 
peak of the pair correlation function. We restrict our 
interest to the monodisperse system where each particle 
has an identical diameter ctq. 

Let us discuss the spatial correlation of the density 
field. Figures [TH] and [1^] present the isotropic parts of 
the structure factor S{k) and the pair correlation func- 
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FIG. 17: (Color online) Collapsed data of the shear rate de- 
pendence of the shear stress S in the nearly elastic monodis- 
perse system (A^ — 2000) with A = 1 using the scaling law, 
Eq. ((5]), for D = 3. The the dotted line and the solid line are 
proportional to 7^ and 7^^, respectively. The legends show 
the volume fraction </!>. The critical exponents estimated from 
Eq. (fT5)l are y^, — 1 and = 2/5. 



tion g{r) for the three-dimensional monodisperse system, 
respectively. Here, g{r) and S{k) respectively satisfy [3l| 



air) 



^EE'^(---^.))' (22) 



S{k) = 1 + (2^)^/2^ 



dr r^-i 



, ^ JD/2-i(fcr) 
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FIG. 18: Structure factor S{k) in the three- 

dimensional monodisperse system {N — 20000) with 
(t> = 0.629, 0.639, 0.649 and 7 = 5 x 10^" and 5 x 10"^ 



g(r) 








■4- 








5x10 






= 0.629 






4 




- - - Y = 


5x10 




* 


= 0.639 






4 






5x10 




<t> 


= 0.649 




5x10 


5 


•f 


= 0.629 






5 






5x10 






= 0.639 




5x10 


5 




= 0.649 



3 

r 



where JD/2-i{kr) is the Bessel function, Sd is the surface 
area of the D-dimensional unit sphere given by Sd — 
27r^/Vr(D/2) with the Gamma function r(£i/2). We 
use A = 1, the viscous contact force given by Eq. ([2]), 
and the number of the particles N = 20000. Because the 
isotropic parts of the spatial correlations are dominant in 
our system, we only focus on the critical properties of the 
isotropic parts in this section. The first peak of g{r) is 
larger than 3 which is the maximum value of the vertical 
axis of Fig. [TH From Figs. [18] and [TH it seems that 
there is no obvious dependence of S{k) and g{r) on 7 
and (j), but the height of the first peak go of g{r) strongly 
depends on 7 and cf). We plot the dependence of the first 
peak of g{r) on 7 for small r — ao region in Fig. 1201 in 
which the first peak go becomes higher and the width of 
the half-height of the first peak ho becomes narrower as 
7 decreases. 

The dependence of the peak on </> and 7 can be roughly 
estimated from our scaling law for the pressure P in Eq. 
([6]). Indeed, the coordination number Z and the pressure 
P, respectively, satisfy 

Z = ^ £" drr^-'g(r), (24) 

p ^ drr^{ao-r)^g{r). (25) 

The derivation of these equations is briefly explained in 
Appendix B. Since the first peak near ao is characterized 



FIG. 19: Pair correlation function g{r) in the three- 
dimensional monodisperse system (TV = 20000) with (j) ~ 
0.629, 0.639, 0.649 and 7 = 5 x 10"'' and 5 x 10"^ 



by the peak value go and the width ho, Eqs. ([M]) and 
(I25p are approximated by 

Z San / r^^^go 

J (jQ-ha 

^ goho{l + 0{ho)}, (26) 
P - Son'' / drr^fc(ao - r)^(?o, 

J (Jo -ha 

^ h^+^go{l + 0{ho)}. (27) 

With the aid of Eq. (pS)) . we find that go and ho satisfy 
the relation near the point J 

50^0 ~ const., (28) 

where we have used the known result on the coordination 
number Z ~ 2D for the frictionless granular particles 
near the point J. Substituting Eq. (^5)1 into Eq. (P7)) . we 
obtain the relation between 170 and the pressure P 

go ^ P-^/^. (29) 

From Eqs. (HH), USD and the height of the first 

peak value go in the unjammed phase is given by 

50^7-2/A|$|4/A^ (30) 
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FIG. 20: The first peak of pair correiation function g{r) in 
tlie tliree-dimensional monodisperse system (A^ = 20000) with 
A = 1 and (j> — 0.639 for various shear rates 7. 



Similarly, from Eqs. (10), p3)) and ((29)) go satisfies 

90 |$r' (31) 

in the jammed phase, which is consistent with the nu- 
merical result of the unsheared jammed system [l^. At 
the critical point, i.e. $ = 0, thus, go is given by 



50 ~ 7 



(32) 

from Eqs. (11), and ((29)) . 

In order to check our predictions in Eqs. ((30)) -((32)). 
we plot go versus 7 for various densities in Fig. [2T] 
The numerical data are consistent with the theoretical 
prediction in which go is proportional to in the 

unjammed phase (see Eg. ((30)) ). but go is almost a con- 
stant in the jammed phase (see Eq. ((3T)) ). and go satisfies 
c,o ~ 7-2/(A+4) pqJjj^ j (ggg Eq.((32])). 

Figure [22l examines the validity of Eq. ((3T)) in the zero 
shear limit of the jammed phase, in which 170 seems to 
satisfy go ^ 1/1 'I' I- On the other hand, go tends to satisfy 
3o7^^^ |<1>|^/* in the unjammed phase as in Eq. ((30)) 
(see Fig. 



FIG. 21: The height of the first peak go of g{r) as a function 
of 7 for the three-dimensional monodisperse systems (A'^ = 
20000) with A = 1 for various densities <j!>. The critical density 
is <pj = 0.639. The sohd line is proportional to 7-2/(^^+4). 
The dashed line is proportional to 7"^''^. 
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FIG. 22: The plots of go versus |<I>| in the jammed phase 
for the three-dimensional monodisperse system (A'' — 20000) 
with A = 1 for various shear rates 7. The dashed guide line 
is proportional to l"!"!"^. 



V. DISCUSSION AND CONCLUSION 

This section consists of two parts. In the first part, we 
will discuss our results, and we will conclude our work 
in the second part. In the first part, let us compare our 
results with those by Hatano [2l|, discuss the growing 
length scale near point J, the long-range spatial correla- 
tion, and the scaling laws for the Langevin dynamics and 
the frictional particles. 

A. Discussion 

Let us compare our results with those by Hatano [2l| . 
Hatano estimated the values of the exponents as 



2.5, 
0.63, 



1.3, 
= 1.2, 



y: 



7 



1.2, 
= 0.57, 



(33) 
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FIG. 23: The plots of goj^^^ versus |<J>| in the unjammed 
phase for the monodisperse system (A'^ = 20000) with A = 1 
for various shear rates 7. The dashed guide line is propor- 
tional to I'l'l''/^. 
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for A = 1, and 

Vi = 



3.2, X. 
0.75, 



7 ~ 



1.3, 
= 1.8 



1.8, 
= 0.72, 



(34) 



for A = 3/2. The system analyzed in Ref. [2l| corre- 
sponds to our polydisperse system, but these values are 
different from those in Eq. Here, let us clarify the 

origin of the differences, (i) The estimation of the critical 
exponents strongly depends on the range of 7. The value 
of and the range of 7 in Ref. [2l| are larger than 
ours. If we adopt and the range of 7 in Ref. [^H, 
we can recover his scaling laws in Eqs. ((33)) and p4l) as 
shown in Fig. [24l while we find the obvious violation of 
his scaling (Fig. [25]) in the small 7 region (7 < 10"''). 
On the other hand, our scaling is still valid as shown in 
Fig. [26] for 7 < 10~^. In order to extract the critical 
properties, it is needless to say that we should use the 
data in the small 7 region. Hence, our prediction for the 
exponents is more appropriate than Hatano's estimation, 
(ii) The estimated exponents in Eqs. (|33p and ((34)) can- 
not be valid even when we adjust the value of as the 
fitting parameter in the small 7 region. The value of 4>j is 
well established by the simulation of the sphere packing 
for the bidisperse system and the monodisperse system 
[isL [l9l | . The scaling plots using the estimated value of 
(Figs. OandO evidently support the validity of our 
prediction. 



{yif 100 



: 0.500 ° <l> = 0-650 

: 0.600 » <l> = 0.660 
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■■ 0.645 



FIG. 24: The scaling of the shear stress S/\^\^'* for the 
polydisperse system {N = 4000) with A = 1 using Eq. p3p 
as a function of the scaled shear rate 'y/|<I>|^'#'/^^ for 10"'' < 
7 < W~\ 0.5 < <;/> < 0.7 with D = 3. 



There are some studies to indicate the diverging time 
scale near point J [H, [2^. Indeed, in the scaling rela- 
tions of T, S, P, and cj in Eqs. ©-((T]) with Eq. the 
shear rate 7 is scaled by the same diverging time scale 
T ~ |c()|-('^+4)/2_ Pqj- conventional critical phenomena, 
the divergence of the time scale is associated with the 
diverging length scale. Therefore, there are many papers 
to discuss growing length scale in the vicinity of the jam- 
ming transition 0, H M, 0, [O, [H, [H H HI . However, 
the length scale in the previous studies is as much as the 
size of several particles in their papers, and it is not clear 
whether the length scale diverges at point J. We, thus. 
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FIG. 25: The scaling of the shear stress S'/|<E>|''* for the 
polydisperse system (A'^ — 4000) with A = 1 using Eq. (|33p as 
a function of the scaled shear rate 'y/l^lf -^/f^ for 5 x 10"^ < 
7 < 5 X 10"^ 0.624 < < 0.652 with D = S. 
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FIG. 26: The scaling of the shear stress S/l^j^* for the 
polydisperse system (A*' — 4000) with A = 1 using Eq. (|15p as 
a function of the scaled shear rate 7/|<I>|^*/!'t for 5 x 10^^ < 
7 < 5 X 10-\ 0.624 <(p< 0.652 with _D = 3. 



conjecture that the length scale might be unrelated to de- 
termine the critical exponents. The success of our mean 
field theory supports the validity of this conjecture. 

In relatively dilute sheared systerns, we know the exis- 
tence of the long-range correlation [H, [13, [H, HI] . The 
existence of the long-range correlation is only confirmed 
in the relatively dilute systems such as (j) < 0.50 for 
the three dimensional case. We, thus, still do not know 
whether there is the long-range correlation in the jammed 
systems, and the role of the correlation. We will discuss 
the long-range correlation in the jammed systems else- 
where. 

The similar scaling relations with different values of the 
exponents are observed in the zero temperature limit of 
Langevin thermostat system, where the Newtonian law 
5 cx 7 is held in the unjammed region |6||. However, we 
cannot extend our simple theory to this system because 
the characteristic frequency defined by a; = 'jS/{nT) is 
always constant in this system, which differs from the 
scaling relation ([7]). We will discuss the results of this 
situation elsewhere. 
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In this paper, we restrict our interest to the frictionless 
particles. When the particles have friction, the situation 
will be changed completely. For examples, the critical 
density (pj of the frictional particles in the static granu- 
lar packings, becomes smaller than that of the friction- 
less particles [37j . and depends on the packing process 
[2^ . [3^. Thus, critical properties of the scaling in the 
sheared dynamical systems of frictional particles should 
differ from the frictionless assemblies. This also will be 
our future work. 



B. Conclusion 

In conclusion, we have extensively checked the validity 
of the mean field theory proposed in Ref. [1^ numeri- 
cally, and we demonstrate that most of all our numerical 
results are consistent with the theoretical predictions in 
Eq. (fT5|) . Thus, we may conclude that the jamming 
transition for frictionless sheared granular materials is 
a continuous transition in which the critical exponents 
are independent of the spatial dimension and are deter- 
mined by the local elastic force between contacted grains. 
We also confirm that the viscosity diverges at point J as 
{(pj — Essential new findings beyond Ref. [l^l are 

the critical behaviors of the first peak of the pair corre- 
lation function given by Eqs. (PTjl and which 
are also consistent with our numerical simulation. 
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APPENDIX A: DERIVATION OF THE VALUES 
OF THE CRITICAL EXPONENTS 

In this appendix, we theoretically determine the val- 
ues of the critical exponents in Eq (|15p . The derivation 
is parallel to that in Ref. [1^, but contains some gener- 
alizations with the help of a simplified argument. 

At first, we should note that the scaling functions 
T±{x), S±{x), V±{x), and yV±{x) satisfy 



lim T_(x) 



X , lim S- (x) 



lim T+(x) 



lim S+{x) = 1, 
2:— *o 



limP_(x) = x^, lim >V_(a-) = x, (A2) 



lim T±{x) ^ x""^, lim S±{x) = x^^, 

lim -P±(.T) = x^'-', lim W±(x) = a;'^^. (A3) 

X — >oo X — 'OO 

The scaling relations P^ - P^ are obtained from Eqs. 
(IXll)-(E3) with Eqs. O-dZi). 

Let us assume that the inverse of the shear rate 7"^ 
can be scaled by a characteristic time scale. Therefore we 
can assume that the ratios between the exponents x^/x-y, 
yct>/y7, y'rj,/y'f and z^/z^ in the scaling laws Eqs. d!])-© 
are common as 



X0 
x. 



7 



y'. 



Z£ 



(A4) 



In other words, the characteristic time scale t exhibits 
critical slowing down as r ~ |<I'|~". This property has 
already been indicated by Ref. [2l|. It should be noted 
that we can determine the exponents without this as- 
sumption (23j . 

Substituting Eq. ^ into Eq. ^ with Eq. ^K^, we 
obtain the relation between the exponents as 



z^ ^ y^ - x^ + a. 



(A5) 



Let us assume that the pressure in the jammed phase 
(<& > and 7^0) converges to that of unsheared 
jammed phase satisfying P ~ 19]. Hence, comparing 
the scaling property of P in Eq. '^'^ T,,;fi^ p iT.a 
find the relation 



(Uni) with P - we 



J/; = A. (A6) 

We also assume that Coulomb's frictional law is held 
in g ranular systems ^2^- Thus, S/P is independent of $ 
[23| and we obtain 

y<i> = y'.p^ (A7) 

with the aid of Eq. p3)) . 

We can use the properties of the cut off frequency /c in 
the density of state in the jammed phase, which satisfies 
fc VP [2(/] . Since we expect that there is only one time 
scale, it is reasonable to assume that the characteristic 
frequency to in the limit 7 — + can be scaled by the cutoff 
frequency fc- Thus, we obtain 



(AS) 



Finally, let us use a similar argument on the character- 
istic frequency lo in the unjammed phase to that in the 
jammed phase. Since the characteristic frequency lo is 
estimated as oj '--^ ^jT/m/l{^) with the mean free path 
Z($), which is evaluated as (cr/Z?0j)|$| in the vicinity of 
point J, we obtain 



\miV+{x) = 1, lim>V+(x) = l, (Al) 

x— )-0 a:— >0 



Xtf, = 2z0 -I- 2, 



(A9) 
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where we have used Eq. ((12)) . 

From Eqs. (|A5p - (|A9p . the exponents a, x^, y^, y'^, 
and Z0 are given by 



= A, = y 



x$ = 2 + A, j/$ = A, 
A 



(AlO) 



The exponents x-y, y^^, y' and obtained from Eqs. 

(jA4p and (jAlOp . Hence, we have determined aU the crit- 
ical exponents as Eq. (fT5|) . 



APPENDIX B: THE EXPRESSIONS OF P AND Z 

BY y(r) 



In this appendix, we derive Eqs. (|24|) and (I25p . The 
coordination number Z is given by Z = M/N , where M 
is the number of the points of contact. Since M is given 
by the number of the pairs whose distances are smaUer 
than (To, Z is given by 



Substituting Eq. (|22p into this equation, we obtain Eq. 



Since the contribution to the pressure P from the elas- 
tic force is dominant in our system, we approximate P in 
Eq. (dni) with the aid of Eq. dUD as 



2DV 
Son 



drr^ fci{r)g{r). 



(B2) 



Z = — I dr l-^'^5{r - r^j)\ . (Bl) Substituting Eq. ([T]) into Eq. (|B2]), we obtain Eq. ^ 
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